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ABSTRACT 



A finite element method was applied to the unsteady- 
transonic small disturbance equation and integrated until 
the solution converged to the steady state for a thin non- 
lifting airfoil. The method of weighted residuals was used 
to formulate the finite element equations, and Houbolt's 
method of central differencing in time was used to integrate 
these equations. 

A secondary investigation applied the steady transonic 
small disturbance equations to a converging-diverging nozzle. 
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I. INTRODUCTION 



Transonic inviscid flows past a smooth airfoil may be 
expressed in terms of the velocity potential ‘i’ satisfying 
the transonic small disturbance equation, 

(1 - mI - M^Cy * * '•’yy “ “ 

This equation presents two major difficulties, 1) it is non- 
linear and 2) it is of the mixed hyperbolic- elliptic type. 
Analytical solutions to non-linear equations are difficult 

to obtain. One must normally resort to numerical methods. 

2 2 

When the coefficient (1 - (y + equation 1 

is negative, the flow is supersonic and the equation is 
called hyperbolic; otherwise the flow is subsonic and the 
equation is elliptic. The forms of the two solutions are 
fundamentally different. Hyperbolic equations admit both 
discontinuities, which propagate only in characteristic 
directions, and the presence of shock fronts. Elliptic 
equations, on the other hand, require that the dependent 
variables and their derivatives be continuous and that a 
change in any part of the flow field affects every other 
part. Many non-linear elliptic equations are solved by 
appropriate relaxation iteration techniques by casting the 
equation in Poisson's form with the non-linearity acting as 
the driving force. Solutions to hyperbolic equations are 
usually obtained by the method of characteristics or by 
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finite difference marching techniques which use an artificial 
viscosity to represent the average jump conditions across the 
shock wave. In mixed supersonic and subsonic flows, normal 
numerical procedures break down because the boundary between 
the two regions is not known a priori. 

Finite element numerical techniques have evolved as a 
powerful tool in obtaining approximate solutions to a wide 
variety of engineering problems, particularly ones with 
Neumann-Dirichlet boundary conditions, i.e., elliptical prob- 
lems. They offer several outstanding advantages. Some of 
these are: 

1) Non-homogeneous problems may be treated with relative 
s implicity . 

2) Complex geometries may be modeled with relative ease 
since the elements can be graded in size and shape to 
follow boundaries of arbitrary shape. 

3) Once the finite element model has been established, 
a variety of problems can be solved by supplying the 
computer with the appropriate data. 

Chan and Brashears [Ref. 5] developed a finite element 
computer program for steady transonic flow over a non- lifting 
airfoil. This program uses the least squares method of 
weighted residuals to approximate equation 1 by a system of 
algebraic equations, and assembles the equations in a special 
way to account for the hyperbolic region of flow. This tech- 
nique prevents the influence of downwind nodes from propagating 
upstream in the supersonic region. 

The purpose of this thesis is to investigate the possi- 
bility. of speeding the convergence of a solution by transform- 
ing the steady transonic equation to the unsteady equation 
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and integrating over time until the time dependent terms 
vanish and to extend the program of Ref. [5] to the tran- 
sonic region of a converging-diverging nozzle. 



II. DISCUSSION OF THE FINITE ELEMENT APPROACH 



In a continuum problem of any dimension, the field 
variable, whether it is velocity potential, velocity, temper- 
ature, displacement or some other quantity, possesses infin- 
itely many values because it is a function of each generic 
point in the solution region. Consequently, the problem is 
one of an infinite number of unknowns. The finite element 
approach subdivides the solution domain into a finite number 
of subdomains called elements and expresses the unknown field 
variable in terms of assumed approximating functions within 
each element. The approximating functions are sometimes 
called interpolating functions and are defined in terms of 
the values of the field variables at specified points called 
nodes or nodal points. Nodes usually lie on the element 
boundaries where adjacent elements are considered to be con- 
nected. In addition to boundary nodes, an element may also 
contain interior nodes. The nodal values of the field variable 
and the interpolating function for the elements completely 
define the bahavior of the field variable within the elements. 
Once these unknowns are found, the interpolating functions 
define the field variable throughout the assemblage of the 
element . 

Clearly, the nature of the solution and the degree of 
accuracy of the approximation depends not only on the number 
and size of the elements used but also on the interpolating 
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functions selected. Interpolating functions may not be 
chosen arbitrarily because certain compatibility conditions 
must be satisfied. Often such functions are selected so 
that the field variable or its derivatives are continuous 
across adjoining element boundaries. Once the problem is. 
formulated in terms of individual elements, the contributions 
of each, element may be assembled to define the entire solu- 
tion domain. This means, for example, that if we are treating 
a problem in stress analysis we can find the force- displace- 
ment or stiffness characteristics of each element and then 
assemble the elements to determine the stiffness of the whole 
structure. Finite element solutions are not, of course, re- 
stricted to structures problems, but the matrix of equations 
defined by the interpolating functions and the nodal field 
variables is still referred to as the stiffness matrix 
regardless of the field variable in the problem. 

Solutions to continuous problems by the finite element 
approach always follow a systematic step-by-step process. 

This process is completely general to the finite element 
method and it is outlined below. [Ref. 4] 

1. Discretize the continuum. 

The first step is to divide the solution domain into 
elements. A variety of element shapes may be used and one 
or more different element shapes may be used in the same 
region. The type and number of the elements used in a given 
problem are a matter of engineering judgement. 
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2. Select the interpolating functions. 



The next step is to choose the type of interpolating 
function to represent the variation of the field variable 
over the element. The field variable may be a scalar, a 
vector, or a higher order tensor. Often polynomials are 
selected as interpolating functions for they are easy to 
integrate and differentiate. The degree of polynomial 
chosen depends on the number of nodes and the nature and 
number of the unknowns and the continuity requirements 
imposed at the nodes and the element boundaries. The un- 
known quantities at the nodes may be assigned to the field 
variable and their derivatives. 

3. Find the element properties. 

After the elements and. their interpolating functions 
have been selected, the matrix of algebraic equations which 
express the properties of the individual elements must be 
determined. Several methods are available for this task. 
These are; the direct approach, the variational approach, 
the method of weighted residuals and the energy balance 
approach. Reference [4] is a good source of information 
on the various techniques. 

4. Assemble the element properties to obtain the system 
equations . 

To find the properties of the over-all system, the matrix 
equations expressing the behavior of each element must be 
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added to the matrix equation of all other elements. The 
basis for this assembly procedure stems from the fact that 
connecting elements have common nodes and the field variable 
must be the same for each element sharing that node. 

At this point the boundary conditions for the system of 
equations must be accounted for and the system of equations 
must be modified before it is ready for solution. 

5. Solve the system of equations. 

The assembly process of step 4 produces a set of simul- 
taneous equations which can be solved to obtain the unknown 
field variables. Linear equations have a number of standard 
solution techniques readily available, but solutions to non- 
linear equations are more difficult to obtain. 

6. Make additional computations if desired. 

Important parameters, such as pressure coefficient in 
aerodynamic problems, may now be calculated from the values 
of the field variables. 
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III. THEORY AND BASIC EQUATIONS 



A. STEADY TRANSONIC FLOW 

Chan et al . [Ref. 5] developed an algorithm to analyze 
steady transonic flow over non-lifting thin airfoils. Bound- 
ary layer effects were ignored and the imbedded shock wave 
was assumed to be weak. These assumptions are consistent 
with transonic small disturbance theory which can be expressed 
mathematically by the following expressions. 

Cl - - m^Cy . * 'fyy = “ CD 

Boundary conditions - 

V • (j) = 0 at infinity (2) 

V = (1 + u)dg/dx on the airfoil (3) 

’u = 0 on the line of symmetry 

( 4 ) 

where g(x,y) defines the airfoil and dg/dx describes the 
airfoil slope. 

The above expressions appear in their dimensionless form 
where <|> = perturbed velocity potential and the perturbed 
velocity components in the x and y directions are respectively 
defined as 

u = 

X 

V = 

y 
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= freestream Mach number and y = ratio of specific heats 
which for air is taken to be 1.4. The physical coordinates 
x' and y' and the velocity potential <p ' are related to the 
dimensionless quantities by 



X = x'/c, y = y'/c, = (})'/cU^ 



where c is the chordlength of the airfoil and U is the free- 

oo 

stream velocity. 

Once the flowfield solution is determined in terms of the 
perturbed velocities, the secondary unknowns are subsequently 
calculated. These include: 



= rl-Im2 



U 



OO 2 1 -i- 



^ = tV(u: - * Csr) ] 






■^0 



[l * 



y/ (y + 1) 



'Y = 



2u . .,2-« u^ , v^ 

JT " -I " -Z 

“ u_ u„ 



2u 

U 



(5) 

(6) 

(7) 

( 8 ) 



In the above, = 1 is the normalized freestream velocity, 

a = local sound speed, p = local static pressure, M = local 

Mach number, V = total velocity, Pq = stagnation pressure 

and C = pressure coefficient. 

P ^ 

B. UNSTEADY TRAiNSONIC FLOW 

Unsteady transonic inviscid flow may be expressed in terms 
of the velocity potential <P(x,y,t) to a first order approxima- 
tion by 
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0 



(9) 



(1 - (y + 1) ^ ) <t> + <|) - = 

'-- oo oo''' XX yy t tt 

which has the same non-linear coefficient retained for steady 
transonic flow in Equation 1. 

Boundary conditions require that the disturbances vanish 
far from the airfoil, 

4> = 0 = 0 at infinity (10) 

X y 

and that the flow remain attached to the body. Let B(x,y,t) = 0 
define the body at any instant. The surface tangency restraint 
may now be expressed by the substantial derivitive DB/DT 
vanishing. 

DB/DT = B^ + (1 + (|)^)B^ + 4>^By (11) 

If the body remains stationary B^ = 0, and the tangency 
condition becomes the same as in steady flow, namely 

V = (1 + u)dg/dx (12) 

where dg/dx represents the airfoil slope. 

The pressure coefficient for isentropic unsteady com- 
pressible flow is defined by 

yM 

Expanding by the binomial expansion and retaining only the 
first order terms gives, 

Cp = -2$^ - 2*^ (145 



{ 



[1 



izl 






(2 + 2 



,^2)]Y/(Y-1).i- 

(13; 
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iAllf. 




IV. FINITE ELEMENT FORMULATION 



The method of weighted residuals is a technique for 
approximating solutions to linear or non-linear partial dif- 
ferential equations and it is the basis for the finite element 
formulation of the transonic small disturbance equation 
(Equation 1) . This procedure involves assuming the general 
functional behavior of the field variable which would approx- 
imately satisfy the basic equation and boundary conditions. 
Substituting this approximation into the original differential 
equation results in some error called a residual. A system 
of algebraic equations results when a weighted average of the 
residual is forced to vanish as it is averaged over the entire 
domain. 



A. STEADY FLOW 



The approximate solution to equation 1 is assumed to be 

$ = N.(j). (15) 

where .N^ are the interpolating functions which exhibits the 
behavior of equation 1 and (j)^ are the undetermined parameters 
at the nodal points. 

When $ is substituted into equation 1, the resulting 
residual is 



R = [1 - M' 

*• o 



m^Cy 

rrk • 



+ 1)N, d), ]N. + N. 

k,x^k-' j,xx j,yy 



C16) 
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The weighted average is determined by multiplying the 
residual R by m linearly independent weighting functions 
and integrating over the elemental domain. Forcing this 
residual to vanish yields, 



Chan et al. [Ref. 5] found that when the weighting function 
for equation 1 is chosen to be 3R/3<j)^ the resulting matrix 
is positive definite and well conditioned. This choice of 
weighting functions is referred to as the method of least 
squares because it is equivalent to minimizing the square of 
the residuals summed over the domain with respect to the un- 
determined parameters. That is, 



ffVf^RdA = 0 



i = 1 to m 



X = //R^dA 



9x/ 3 = //8R/3 <J>^RdA = 0 



(17) 



Integrating over the domain produces the system of 



algebraic equations 



(18) 



where the elemental matrix S^j is defined as 




(19) 



With Qj and equal to 
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B. UNSTEADY FLOW 



Development of the unsteady flow finite element equations 
is similar to the procedure used to formulate the finite ele- 
ment equations for steady transonic flow. The least squares 
method of weighted residuals is again used but <f> is now a 
function of time as well as the spatial coordinates x and y. 
The approximate solution has the form, 

9 = N^(x,y)(j)^(t) (20) 

Substituting ^ in the unsteady transonic small disturbance 
equation, the weighted residual becomes, 

X = //(Rj_ + R 2 + R3)^dA (21) 

where 



R2 = 



R, = 






Expanding equation 21 yields 



X = ffi^i ^ ^ ^ ^ 2R^R2 + 2R2R3]dA (22) 



and minimizing with respect to the undetermined parameters 
(j)^ the following system of algebraic equations result, 

(j)j = 0 = 2//3R^/3(j)j [Rj^ ^2 R3]dA 
OR^/act-j = Pi 
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where has been previously defined in the steady finite 
element formulation 

The above equation may be rewritten in the form 



S . . . + SC . . <j) . + SM ..(}). = 0 

iri iri i; j 



(23) 



The stiffness matrix is the same as that developed 

for the steady transonic equation and the damping (SC^^) and 

mass fSM..) matrices are defined below. 

^ ij 



.SCij = -/ACN.^^P.dA 



SM. . = -//M^N.P.dA 

ij j 1 



(24) 

(25) 
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V. ELEMENT DESCRIPTION AND ASSEMBLY OF EQUATIONS 



A. ELEMENT DESCRIPTION 

The basic element used in the finite element program is 
the non- conforming cubic triangular element developed by 
Bazeley et al . [Ref. 2]. Also used in the program are the 
quadrilateral and trapezoidal elements constructed from these 
triangular elements. These three types of elements can be 
mixed and used freely in the entire flow region except that 
only trapezoids should be used to cover the supersonic and 
mixed region in order to enact the special assembly procedures 
required by the hyperbolic equation which describes the flow 
in that region. 

The basic triangular element is shown in Fig. 1, which at 
each vertex has the velocity potential and the velocity com- 
ponents as the undetermined parameters. This type of element 
was chosen because both Dirichlet and Neumann boundary condi- 
tions can be imposed with equal convenience. In addition, 
because velocity components are used as primary unknowns 
secondary parameters, such as Mach number and pressure coef- 
ficient can be calculated directly without resorting to 
numerical differentiation, which would produce additional 
errors. 

In the element, the approximate solution is assumed as 
(f) = (i ~ 1 P) 
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Figure 1 - Triangular Element 
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In which <j) ^ ' s are the nine undetermined parameters of (p and 
N- are the interpolation functions which are expressed in 
terms of the area coordinates. 

The shape or interpolating functions and their first and 
second derivitives are defined below. 

Letting , 

^i = 

^i = >^j - >^k 
^i = ^^k - 

A = area of triangle 1-2-3 = (b-c, - b,c.)/2 

J K K j 

a = 0. 5 (c^ - c^) 

B = 0.5 Cbj - bj^) 

^ ^ k^ j ^ k 

^ = ^i^j^k ^ ^j^k^i " ViS 
Hxx = * Cjb^b,^ - 

yy ^^ijk ^jki ^kij-^ 

with i = (l,2,3),k = (3,1,2) then one has for 1 = (1,4,7), 
i = (1,2,3) 

Ni = qh3 - 2;.) * 2H 

" “x1/2A 

^l,y " f^^i^i ■ ^i^ ^ 2H^]/2 a 
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f'l.xx - * “xxl 

“^l.yy ' t^'^i ^ 2»yy’ 



C2,5,8), 


, 1 = (1,2,3) 










"’l 


■ 'i f^k'j ■ =; 


j?k) * CH 








"'^l,x 


= [2bi?i(c^?. ■ 


■ '"j k> * 


1 ^} 

1 


+ aH ] 
x-* 


/2A 


i,y 


” [2^i«if^'j ■ 


■ '^j'k' * 


aH ] 
y ■' 


/2 A 




^l,xx 


■ [2»i (Ck'j - 


■^j'k’ * 


4b. (2A)?. + 


aH ] 


N-, 

i,yy 


= [2c/ (c^C. - 


=j'k> * 


“Hyy] 


/(2A) 


2 


(3,6,9), 


, 1 = (1,2,3) 










Nl 


= ;/ (b. - b, 


kS’ * “ 








N-, 

l,x 


' [2bi?i(bj?k • 


■ bk'j[ * 


BHxJ 


/2A 




N-i 

i,y 


= [2c.?. Cb.?^ 


- bk'j) 


* 2C.2 
1 


+ 6H 

y 


] /2A 


Nn 

1 ,xx 


- [2bi [b^^k - 


bk?^J * 


“xxl 


/(2A) 


2 



Nl,yy = [2c/ Cb-t, - b,5 .) * 4 c.( 2 A). * BH^yJ /(2A)2. 
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Quadrilateral and trapezoidal elements as shown in Fig. 2 
are also used in the present program, the former is used in 
the subsonic region and the latter in the mixed and super- 
sonic region. The element matrix for the quadrilateral ele- 
ment is obtained by combining appropriately the matrices for 
two triangles, while the matrix for trapezoidal element is 
obtained by averaging contributions from two left-running and 
two right- running triangles. The averaging process is designed 
to remove the bias effects inherent in the quadrilaterals used. 

B. ASSEMBLY OF EQUATIONS 

Straightforward application of the finite element assembly 
technique to transonic flows would fail (the solution diverges) 
because this would allow disturbances to propagate upwind in 
the supersonic region of flow where the governing equation is 
hyperbolic. Hyperbolic equations have a time-like dependency 
in that the solution at the downwind station is affected by 
the upwind station but not vice-versa. Assembly techniques 
for a transonic flow finite element program must take into 
account this time-like dependency. If the x-axis is taken 
as a time-like direction in the supersonic region, the ele- 
ment matrix may be assembled in a way similar to a backward 
finite difference operator, which has been successful in 
solving hyperbolic equations. 

Consider the rectangular element sketched below with the 
upwind station I and the downwind station II, each having two 
nodal points with the element type chosen. The element 
matrices can be constructed in the usual manner. ' 
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a. Quadrilateral 
Element 




b. Trapezoidal 
Element 



Figure 2 - Quadrilateral and Trapezoidal Elements 
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I 

However, before assembling the element matrix into the system 
matrix the non-linear coefficient of equation 1 is evaluated. 

C = 1 - (y + 1) u 

The sign of the coefficient being positive, zero, or 
negative defines the equation as elliptic, parabolic, or 
hyperbolic. If C is non-positive for all nodes in the element, 
the rows representing the improper - downwind influence on the 
solution at an upwind station are ignored during assembly. 

This feature is automatically applied in the program requiring 
only a little care in arrangement of the nodes of the element. 
In the anticipated supersonic region, element node points 
should be arranged in the order as indicated in figure 3, 
starting with the upper left corner and proceeding in the 
counter-clockwise direction. In the elliptic region, i.e., 
where the coefficient is positive, no special assembly tech- 
nique is invoked. 

C. ITERATIVE PROCEDURES 

With the equations assembled and the proper boundary con- 
ditions imposed, the system of non-linear algebraic equations 
is solved by iterative procedures in the form 

= 1. (23) 

t" Vi 

to solve for the solution in the n^ iteration. The function 
(|) is defined as 
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I 



II 



I 4 



FLOW 

DIRECTION 
> 



Figure 3 - Nodal Arrangement for Supersonic Region 
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in which the under-relaxation factor 0 is in the range 
0<9_<1. For subsonic flow 0 = 1, which is simply a succes- 
sive approximation, yields good results, but it is necessary 
to under-relax somewhat with 0 approximately .5 for super- 
critical flow. Generally, a smaller relaxation factor will 
make the solution more stable but it will tend to slow down 
the rate of convergence. 

Equation 23 is subject to the convergence criterion that 
the change in local Mach number between two successive iter- 
ations is less than a prescribed value e at all nodes in the 
flow field. That is, 



,,Cn) . ^(n-1) 



< e • 



/ 
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VI . INTEGRATION OF UNSTEADY FINITE ELEMENT EQUATIONS 



The unsteady transonic small disturbance equation (equation 
9) when suitably reduced to a finite element approximation 
appears in the form, 



S. .<t> . + SC- .4> . + SM. .(}) . = R. 
II 3 iJ J iJ I 3 



(26) 



This equation is analogous to a damped spring mass system, 
hence S-., SC--, and SM. . are respectively referred to as 
the stiffness, damping, and mass matrices. 

Mathematically, equation 26 represents a system of second 
order differential equations with constant coefficients, which 
can be solved by standard numerical procedures for differential 
equations, such as Runge-Kutta or Milne methods. However, 
this would be a very costly technique if the coefficient 
matrices are very large. In practical finite element analysis 
there are a few effective methods which take advantage of the 
banded matrices usually encountered in finite elements. One 
such method is the direct numerical integration method. 

Direct integration involves a numerical step-by-step 
procedure aimed at satisfying equation 26 only at discrete 
time intervals At apart and not over all time t. Conceptually, 
direct integration is a finite element method in space and a 
finite difference method in time. Examples of direct inte- 
gration are the central difference method, Houbolt integration 
and the Wilson method. The first two schemes are finite 
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difference schemes whereas the latter is a linear accelera- 
tion method. Linear acceleration integration assumes a 
linear variation of acceleration from time t to time t + At. 

Central differencing can be very effective in the solu- 
tion of many dynamic problems especially those that involve 
a large system of equations. However, this method is unstable 
for all time steps larger than a critical time step. 

Houbolt integration is an implicit finite differencing 
method related to central differencing, only it has the ad- 
vantage of being stable for all TIME STEPS. 

The Houbolt method was used to integrate the unsteady 
finite element equations because of this stability. 

Houbolt integration uses the following finite difference 
expansions : 



n,t+At ' ^*1,1 * “♦l.t-At' ♦i,t-2At]l/'^’'' 












which are two backward-difference formulas with errors of 
2 

order (At) . 

The solution of ^_j_^^must satisfy equation 26 and at 
time t+At equation 26 becomes 



S . . d) . , SC . . i . , SM . . d) . . ~ 0 

ij^;),t+At i;)^j,t+At ij^j,t+At 



Substituting the finite difference formulas for (p . and 

rearranging all known vectors on the right hand side, the 

solution for <p . _^,^is obtained, namely: 
j , t+At ’ 
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Ca^SM. . 

0 ij 



a, SC. . 
1 ij 



+ 



s . O'f’ . 

ir j,t+At 



, t+At 



(27) 



^ Ca2SM.. + a3SC. + (a4SM.. + a5SC..)<t,. 

Where the constant integration coefficients are: 
a^ = 2/At^ 

= ll/6At 
E2 = 5/At 
3. ^ / A ^ 

*4 = 2a„ 

^5 ° ‘^ 3/2 
^6 ■ 

a? = 3 j/9 

Equation 27 may be written as 



, t- 2 a t 



SE. . (p. = RE. 

13 3 > t+At j 



(28) 



where the effective stiffness matrix SE . . and the effective 

13 

load vector RE. are defined as: 

3 / 

SE . . = S . . + a„SM. . + a, SC . . 

13 13 0 13 1 13 

RE. = SM-.Ca^t.,- ^ * 3.. <p. + a.P. . -_) 

3 13^ 2 ^ 3 , t 4 ' 3 ,t-At 6 3 ,t- 2 At-' 

+ SC..(a.,4>. * + ac'1>-. ^. + a_4>. * ^ *) 

13^ 3 »t 5 3.t-At 7 3 ,t- 2 At-" 



are 



Accurate knowledge of the vectors d 

j,t-At 

required to yield an accurate solution for 



and 

3 



P 

3 , t- 2At 



t + At 



and 
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normally the Houbolt integration scheme requires a special 
starting procedure to determine the initial two vectors 
At 2At* However, since the primary interest of 

this problem is to integrate the equations until they con- 
verge to a steady state solution, it is not necessary to 
obtain an accurate time history of the flow. Errors induced 
by the inaccurate starting vectors will vanish as time 
approaches infinity. Therefore, the starting vectors may be 
chosen somewhat arbitrarily. 
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VII. CONVERGING-DIVERGING NOZZLE 



S. F. Shen [Ref. 10] demonstrated the feasibility of 
calculating compressible flows through a converging-diverging 
(Laval) nozzle by dividing the region of calculations into 
three patches, a subsonic region, a supersonic region and a 
transonic one, of course bounded by the other two regions. 

The locations of the boundaries for each region were chosen 
arbitrarily provided the sonic line is bracketed by the 
subsonic and supersonic boundaries. 

Two different finite element formulations were used for 
the subsonic and supersonic regions, but Shen [10] resorted 
to analytical approximations to cover the transonic patch. 

This restricted the calculations to nozzles with small throat 
curvatures because no analytic solutions exist for nozzles 
with large throat curvatures. It is conceivable that STRANL-II 
could be adapted to provide a continuous solution throughout 
all three regions. 

Outside the transonic region of flow the governing small 
disturbance equation is 

(1 - M^)(J> +4) =0 (29) 

V. oo-'^xx ^yy 

This holds for both subsonic and supersonic flow. Comparing 
equation 29 with equation 1, the transonic small disturbance 
equation, we notice that only the non-linear coefficient 
(Y + l)u distinguishes the two equations from each other. 
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This coefficient becomes negligible when the Mach number 
becomes less than .8 or greater than 1.2. With this con- 
sideration in mind, it was assumed that equation 1 would 
adequately describe the flow through the nozzle and that 
the finite element formulation developed for the non- lifting 
airfoil would apply to the Laval nozzle. 

A. BOUNDARY CONDITIONS 

Two solutions are possible for a converging -diverging 
nozzle: 1) Symmetric flow, where the flow is subsonic through 

the domain, except for a small supersonic region near the wall 
in the throat, and 2) Asymmetric solution, where the flow 
accelerates to sonic velocity in the throat and then continues 
to accelerate to supersonic velocity in the diverging section. 
Different boundary conditions apply for the two solutions. 

For the symmetric case, both inlet and exit velocities must 
be specified. Inlet and exit velocities are equivalent in 
the subsonic solution. The supersonic solution requires that 
only the inlet velocities be specified. If the exit veloci- 
ties are also applied, the problem is overspecified and the 

/ 

solution may not ' converge . 

Velocities at the inlet and exit are not uniform in the 
y direction, therefore the disturbances cannot be set to zero 
as in the case of the non- lifting airfoil. Boundary velocities 
must be calculated by solving equation 29 analytically. 

Equation 29 is a linear equation which can be mapped to 
Laplace’s equation, 

(f) = 0, 
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by letting y' (1 - M^) y. Laplace’s equation is easily 

solved for the case of the hyperbolic nozzle by transforming 
from cartesian coordinates to elliptic coordinates. This 
transformation simplifies the solution because the stream 
lines must be hyperbolas to follow the nozzle boundary and 
therefore follow the hyperbolic coordinate v = constant. 

If the elliptic coordinates y and v are chosen such- that 
the curves y = constant are ellipses and the v curves are 
hyperbolas, then the velocity potential which satisfies 
Laplace’s equation for a hyperbolic nozzle is simply 

({) = Ay 

where A. is a constant of integration. The stream function is 

= Av 

The transformation w = y + iv = cosh ^(2z/a) gives rise 
to the elliptic coordinates 

y = 1/2 a cosh y cos v, x = 1/2 a sinh y sin v 

r = 1/2 a [coshy - sin“v] 

= aJ{j + a/2)^ + x^ 

= 'sjiy - ^ y} 

Solving for y and v produces 

y = cosh'^ [ (r^^ + r ^ ! o.\ 

V = cos'^ [Ctj_ + r2)/a] 

The nozzle boundary is defined by = constant, which 
along with the equation for the nozzle wall in cartesian 
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coordinates, y' = £C^), implicitly defines the constant a. 
Substituting for y in the velocity potential produces, 

= A cosh ^ t n2)/a] 

from which the velocities may be determined. 




The constant of integration A may be determined by speci- 
fying the flow rate through the nozzle, but when the inlet 
velocities are normalized with respect to the freestream 
velocity (U^) A is factored out of the problem. 

Velocities for compressible flow can be solved by mapping 
back to the physical coordinate system (x,y plane) . 

Other boundary conditions are universal to both problems. 
These are: 

on the line of symmetry 



u = 0 

V = (1 + u) df /dx 
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at the nozzle wall 



F Cx) defines the nozzle boundary in terms of a ratio of the 
throat semi-height as a function of x. The throat semi-height 
is taken to be 1 for convenience. 

Pressure ratio, sound speed, and Mach number are calculated 
as before by equations 5 through 8. 
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VIII. DISCUSSION OF RESULTS 



A. TIME INTEGRATION TO A STEADY STATE SOLUTION 

As stated before, the Houbolt method of integration is 
stable for all time steps. Results of the test cases bear 
this out with the larger time steps providing the most rapid 
convergence to a steady state solution. Time steps were 
tried from t = . 1 to t = 100. Time steps larger than t = 100 
were not attempted because as t becomes too large the influ- 
ence that the damping and mass matrices have on the effective 
stiffness matrix becomes negligible compared to the stiffness 
matrix. That is: 

SE.. = S.. + 2SC../At^ + 6SM../At 
ij ij 11 11 

as At ->■ 00 



The starting solutions were chosen somewhat arbitrarily. 

d) . - ^ was chosen to satisfy the first iteration of the steady 

^l,2At ^ 

solution 



S 



j , 2At 



0 



when the non-linear term (u) in the coefficient 
1 - - M^(y + 1) u 

was set to zero. *• . and ^ were chosen as multiples of 

1 > At ,0 

Ip ■ ^ ^ and respectively they were 
l,2At 
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This starting procedure proved to be superior to chosing 
the first three vectors closer to the converged solution. 

2At’ At’ 0 chosen to be the last three 

time steps of the previous case, the solution oscillated and 
converged much slower than with the starting solutions chosen 
as above. 

The stiffness, mass, and damping matrices were recalculated 
after each time step, using the under- relaxation technique 
described above. This was necessary to utilize the special 
assembly procedures invoked by STRANL-II to prevent the in- 
admissahle influence of downwind nodes from propagating up- 
stream in the supersonic region. 

For barely critical flow (M^= .861) and subsonic flow, 
an under-relaxation factor 9=1 (successive approximation) 
resulted in convergence to a steady state solution after only 
three time steps. Eleven time steps were required for the 
supercritical solution to converge using the same relaxation 
factor. Reducing 9 to . 5 increased the rate of convergence 
and the solution achieved steady state after six time steps. 
Figure 4 compares the steady state solution for a 61 thick 
circular arc airfoil at M'^ = .909, using the same integration 
method, with the results obtained in Ref. [5]. Chan's results 
converged in 10 iterations after using the results from the 
barely critical flow as an initial guess to the 



supercritical solution. Figure 4 is a plot of local Mach 
numbers at boundary nodes on the airfoil. 

B. CONVERGING-DIVERGING NOZZLE 

The nozzle chosen for the test cases was the two-dimen- 
sional Oswatitsch nozzle with the boundary defined as 

y ~ ^ ^ 1 (x - 2.5)^ 

where the throat at x = 2.5 has a semi-height of 1. The inlet 
was taken to be x = 0 and the exit was at x = 5. M = .44, 
the inlet Mach number on the nozzle center-line was chosen to 
yield sonic conditions in the throat. 

Two solutions were possible for this inlet condition,- the 
symmetric solution and the asymmetric solution; but neither 
solution was achieved by the finite element method. Although 
the solution converged for the subsonic case in three itera- 
tions, center-line Mach numbers deviated significantly from 
both one-dimensional theory and from Oswatitsch' s approxima- 
tion [Fig. 5] . When the local Mach number M exceeded the 
inlet Mach number by approximately .2 (M' ^ .64) the solution 
was invalid. Differences at the center part of the nozzle 
are due to an essentially incorrect free stream Mach number. 
Patching the solution at x - 1.5 would improve the solution. 

A second test case was run for the supersonic section of 
the nozzle with the inlet boundary on the sonic line. The 
exit boundary was left free to float. Here the solution was 
unstable and no meaningful results were obtained. 
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Time integration 
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Figure 4 - Comparison of Time Integration Results with Steady State Results. 



c 
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Figure 5 - Nozzle Center-Line, Inlet Mach Number 0.4300. 



IX. PROGRAM MODIFICATION 



The finite element computer program for non-lifting air- 
voils, as developed in Ref. [5], is separated into two parts. 
These have been designated STRANL-I and STRANL-II by Lockheed 
Corporation. STRANL-I generates a finite element mesh to be 
used as inputs to STRANL-II, which assembles the finite ele- 
ment equations, applies the boundary conditions, and solves 
the non-linear system of equations. Detailed descriptions 
and instructions for the use of the two programs can be found 
in Ref. [5]. Only the modifications to the above programs 
will be discussed in this section. 

A. UNSTEADY EQUATIONS 

Modifications to STRANL-II to form and solve the unsteady 
finite element equations were three-fold: 

1) The new elemental matrices SCj^j and SMj^j were calculated 
and assembled. 

2) All the matrices were stored on an external magnetic 
disk to be accessed and reassembled later because of 
the amount of space required to store three large 
matrices, in core memory. 

3) The effective stiffness matrix SEj^j and the effective 

load vector were assembled, and the system of 

equations solved. 

Several existing subroutines in the original STRANL-II 
program were modified to assemble the damping and mass matrices. 
These include subroutines NEWK, EMTC, DERV, and EMQT. Two new 
subroutines were added to perform the other tasks. 
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B. MODIFICATIONS TO CALCULATE THE MASS AND DAMPING MATRICES 
EMTC in the STRANL-II program calculated the elemental 
stiffness matrix by numerically integrating the equation, 



S.. = Q.P.dA 
11 ^1 1 

Equations were added to EMTC to perform the additional numer- 
ical integrations for the mass and damping matrices. All 
three matrices were calculated at the same time. EMQT 
assembled the elemental stiffness matrices for a quadrilateral 
and trapezoidal element from the contributions of the tri- 
angular elements. Mass and damping matrices were calculated 
in the same fashion. 

Subroutine NEWK, an original subroutine in STRANL-II, 
which assembled the finite, element equations for steady flow, 
was modified to assemble SC. . and SM. . . A new calling argu- 
ment, NMAT was passed to NEWK, which assembled contributions 
from the triangular, quadrilateral and trapezoidal element 
matrices into the global matrices S- - , SC- • , and SM. . , 
depending on NMAT being 1, 2 or 3. 

1 . Subroutine STORE 

Given a non- symmetric matrix stored in a banded node, 
plus the right hand side vector, subroutine STORE separates 
this system into two matrices and stores them on a magnetic 
disk. Figures 6 and 7 show the decomposition of a banded 
matrix into banded storage, and the further decomposition 
of this banded stored matrix to two smaller matrices by sub- 
routine STORE. In these figures, D, L, and U represent the 
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NxN NX2HBW 




-f- 
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Figure 6 - Decomposition of a Banded Matrix Plus Right Hand Side 



NX2KBW NXHBW NxHBW 






1 






Figure 7 - Separation of a Banded Stored Matrix by Store 



diagonal matrix, the lower triangular matrix, and the upper 
triangular matrix respectively. HBW is the half bandwidth 
and R is the right hand side vector. 

STORE requires an additional work area one-half the 
size of the originally dimensioned matrix which is to be 
stored. 



2. Subroutine TIME 

Subroutine TIME integrates the system 



S. . + SC. .4) . + SM. A . = R. 
iJ J 11 1 11 1 1 



by Houbolt integration. 

TIME reassembles the three matrices which were stored 

on the magnetic disk to form the effective stiffness matrix 

and the effective load vector. Once this system of equations 

is assembled, a banded equation solver is called to yield the 

solution ford). ^ , . Figure 8 is a schematic flow chart of 

j,t+At 

TIME. In Fig. 8 when L = 1, the lower triangular matrix and 
the diagonal of the effective stiffness matrix are formed by 
adding the appropriate contributions from the stiffness, mass 
and damping matrices. When L = 2, the upper triangular matrix 
is formed in like fashion. 



C. CONVERGING -DIVERGING NOZZLE 

1 . Application of the Boundary Conditions 

Regardless of the type of problem for which a set of 
system equations have been assembled, the equations will have 
the form 
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Figure 8 - Flow Chart of Subroutine TIME 
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STRANL-II 



F igur e 




STR ANL-I I 

Flow Chart of STANL-II Modification 
to Integrate Unsteady Equations 
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in which is annxn matrix and x^^ and are vectors of 

length n. These equations do not take into account the known 
values of x^ on the boundaries. However, for a unique solu- 
tion of the above equation, at least one or more nodal var- 
iables must be specified and must be modified to render 

it non-singular. For each equation i, either x^ or must 
be specified but it is physically impossible to specify both 
x^ and R^ . There are a number of ways to apply the boundary 
conditions to the equations and when they are applied the 
number of equations is reduced. However, it is convenient 
to leave the number of equations unchanged to avoid major 
restructuring of the computer storage. One such method is 
described below. 

If k is the subscript of the prescribed nodal variable, 
the k^^ row and the k^^ column of the original matrix are 

set to zero, is set to 1 and Rj^ is replaced by the known 

value of Xj^. Each of the n-1 remaining terms of R^ is modi- 
fied by subtracting from it the value This procedure 

is repeated for all the boundary values. Of course, when the 
matrix is stored in a banded mode, the algorithm will differ 
from that for the nxn square matrix, but the procedure is 
similar . 

Subroutine BNDRY applies the boundary conditions for 
the modified program. Setting the option parameter I0PTC4), 
in STRANL-II, equal to 1 will call BNDRY which will read the 
boundary velocities and apply them to a banded stored matrix 
by the method described above. 
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When the value of the x-^ is zero, as in the non- 
lifting airfoil problem, the algorithm becomes simpler than 

the above method because there is no need to either set the 
1" h 

column to zero or subtract the value from the 

right hand vector. 
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X. TEST CASES 



In computing the flow field for either the non-lifting 
airfoil or for the Laval nozzle, the following procedures 
were followed: 

1) The desired mesh was sketched with each node assigned 
a number. 

2) Appropriate input cards based on the sketch were pre- 
pared and supplied to STRANL-I to generate the data 
on punched cards for input to STRANL-II. 

3) The above punched cards were supplied to STRANL-II 
with three additional cards as input parameters for 
each case, plus additional cards for the boundary 
velocities, if the nozzle solution is desired. 

4) Results of the finite element calculations are printed 
after each iteration, and the converged solution is 
punched on cards for possible later use. 

A. TIME INTEGRATION TO A STEADY STATE SOLUTION 

Test cases for the integration of the unsteady transonic 
finite element equations were conducted to calculate the 
steady transonic flow over a 6% thick circular arc airfoil. 
These tests were made using the same airfoil, mesh, and free- 
stream Mach numbers as Chan et al. [Ref. 5] published. These 
conditions were chosen to provide a source for comparison of 
the results. 

Freestream Mach, numbers used in these calculations were: 
M^ = .806 (subcritical) 

M_^ = .861 (barely critical) 

M^ = .909 (supercritical). 
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Each case was treated individually with (j) j = 0 used as 
the initial guess for each case, whereas Chan et al. [Ref. 5] 
used zero for the initial guess for = .806 and then used 
the computed results as the initial guess for each subsequent 
case . 



DELT, the value for the time step, is input by a parameter 
specified in columns 41-45 of the second card following the 
title card for each case when the unsteady option (IOPT(6) = 1) 
is selected. 



1 . STRANL- I Program 

a. Input 

Input cards used to generate the finite element 

mesh are listed on the next page. Cards were arranged in 

accordance with Ref. [5], in the following order: 

Title card 
Option card 
Element cards 

Card for the total number of nodes 
Node coordinate cards 

Card for the number of boundary nodes 
Card for the boundary nodes at infinity 
Card for the nodes on the line of symmetry 
Card for the nodes on the airfoil 
Cards for the slope of the airfoil. 

Input cards to STRANL- I for these tests are listed 
on the next tnree pages. 

b. Output 

Output from STRANL- I is in the form of printouts 
and punched cards. Printouts from STRANL- I are listed on the 
following eight pages. 
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2. STRANL-II Program 



a. Input 

Listed on the next three pages are the input cards 
to the STANL-II program. 

b. Output 

The output from this program is in the form of 
printouts for each iteration and punched cards for the con- 
verged solution. Output from STRANL-II for = .909 is listed 
on the following eight pages. 
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B. CONVERGING-DIVERGING NOZZLE 



Two test cases were run for the converging-diverging 
nozzle. The first case was for symmetric flow designed to 
yield sonic conditions in the throat of the nozzle. The 
second case dealt with supersonic flow in the diverging sec- 
tion by starting with the sonic line as the boundary of the 
nozzle mesh. Oswatitisch' s two-dimensional nozzle [Ref. 9], 

Y - 1 2 (x - 2.5)^ with a semi-throat height of 1 at 

X = 2.5 was used in both cases. Boundaries for the subsonic 
nozzle were at x = 0 and x = 5. 

1 . Symmetric Solution 

This problem was analyzed using the mesh shown in 
Fig. 10, which consists of 126 elements and 152 nodes. In 
the second card (the option card) IOPT(4) = 1 indicates that 
non-zero boundary velocities at the inlet and the exit will 
be read and applied by subroutine BNDRY. This option requires 
that the number of inlet and exit boundary nodes be specified 
in columns 36-40 of the next card. Perturbation velocities 
at the boundary nodes follow on the subsequent four cards. 
Subroutine BNDRY reads u and v respectively for the first 
boundary node and then continues reading u and v for each 
inlet and then each exit node in the order specified on the 
appropriate card. 

2 . Supersonic Case- -Diverging Section 

The mesh used for this case is sketched in Fig. 11 and 
input cards to STRANL-II follow on the next page. For this 
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case the options in effect are I0PT(4) = 1 and I0PT(5) = 1 
which cause non- zero boundary velocities to be applied to 
the sonic line only. 
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INPUT TO STPANL-I I 
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Figure 10 - Finite Element Mesh for the Converging-Diverging Nozzle. 
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Figure 11 - Finite Element Mesh for the Supersonic Case 
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